%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                           Extension Program                             %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ---------------------------- Description -------------------------------- 

    % This program calibrate the extended model to the initial balanced 
    % growth path and the ending balanced growth path. The decomposition
    % exercise identifies the contributions of each potential channels. 
    % The program generates Table 10 in the paper and Appendix Table 17-19.

clc
clear all

global eta lambda zeta
global chi_a chi_b rho_a rho_b gamma_a gamma_b mu iota
global phi nu theta
global epsilon beta delta delta_p
global Q alphah alphal mh ml zbar
global ss

options=optimoptions('fsolve','Display','off');
options2=optimoptions('fmincon','OptimalityTolerance',1E-8,'Display','off');

%--------------------------------------------------------------------------
%   Parameters Values from a Priori Information or Direct Estimation      %
%--------------------------------------------------------------------------

eta=0.85*(1/3);       % Capital Share
lambda=0.85*(2/3);    % Labor Share
zeta=1-eta-lambda;    % Profit Share

epsilon=2;            % CRRA Parameter
delta=0.069;          % Depreciation Rate
delta_p=0.94;         % Patent Survival Rate

r=1/(1+0.06);
g=(1+0.0213)^((zeta+lambda)/zeta);
beta=r*g^(epsilon*zeta/(zeta+lambda)); % Discount Factor

Q=[0.872, 1-0.872; 0.017, 1-0.017];    % Type Transition Matrix
ptt=(1/2)*ones(1,2)';
tol=10^-6;

% Find the stationary distribution of idiosyncratic shocks
for j=2:1000
    pptt=Q'*ptt;
    if max(abs(pptt-ptt))<tol
        ptt_star=ptt;
        break
    end
    ptt=pptt;
end

alphah=ptt_star(1); % Share of High Type
alphal=ptt_star(2); % Share of Low Type
ml=exp(-0.355);     % Prod. Ability of High Type
mh=exp(3.196);      % Prod. Ability of Low Type

chi_a=1;            % Applied Research Cost, Scale, Normalization
zbar=1;             % Normalization
nu=0.7;             % Matching Function, Exponent

%--------------------------------------------------------------------------
%                        Setting Data Targets                             %
%--------------------------------------------------------------------------

g_target=(1+0.0213)^((zeta+lambda)/zeta);
rd_sales_h=0.0362; 
rd_sales_l=0.0283;
range_h=11.81;    
range_l=1.924;     
range_avg=alphah*range_h+alphal*range_l;
trans_target=0.2320;
basic_h=0.0420;
basic_l=0.0373;
basic_target=0.0395;

g_target2=(1+0.0222)^((zeta+lambda)/zeta);
rd_sales_h2=0.0315;   
rd_sales_l2=0.0671;   
range_h2=6.307;       
range_l2=1.609;
range_avg2=alphah*range_h2+alphal*range_l2;
trans_target2=0.3704;
basic_h2=0.0461;
basic_l2=0.1146;
basic_target2=0.0611;

%--------------------------------------------------------------------------
%                     Calibrating to the Initial BGP                      %
%--------------------------------------------------------------------------

ss=0.05; % R&D Tax Credit in 1981-1985

% Calibrate key parameters

par_vec=@(vec)par_solver_ext(vec,g_target,rd_sales_h,rd_sales_l,range_h,range_l,trans_target,basic_h,basic_l);
vec_guess=[2 ,   0.00015  , 1  , 5 ,   0.6,  0.2 ,  0.05  ,  0.15 ];

lb=[0,0,0,0,0,0];
ub=[1E6,1E6,1E6,1E6,1E6,1];
[par_opt,V_opt,exitflag]=fmincon(par_vec,vec_guess,[],[],[],[],[],[],[],options2);

if exitflag <=0
    disp('Parameters Not Solved')
end


iota=par_opt(1);    % Management Cost, Elasticity
mu=par_opt(2);      % Management Cost, Scale
gamma_a=par_opt(3); % Step Size of Applied Research
chi_b=par_opt(4);   % Basic Research Cost, Scale
rho_a=par_opt(5);   % 1+rho_a: Applied Research Cost, Elasticity
rho_b=par_opt(6);   % 1+rho_b: Basic Research Cost, Elasticity
phi=par_opt(7);     % Matching Function, Scale
theta=par_opt(8);   % Bargaining Power of the Patent Seller

gamma_b=1.6*gamma_a;

% Compute moments in the initial BGP
guess=[0.9,0.1,0.02,0.001,11,2,1,6,1,2,0.01];
geq_vec=@(input)geq_solver_ext(input,iota,mu,gamma_a,gamma_b,chi_b,rho_a,rho_b,phi,theta);
[answer,value,exitflag] = fsolve(geq_vec, guess, options);
if exitflag ~= 1
    disp('First-Order Optimality is not Achieved')
end

iah=answer(1);
ial=answer(2);
ibh=answer(3);
ibl=answer(4);
omegah=answer(5);
omegal=answer(6);
a=answer(7);
v1h=answer(8);
v1l=answer(9);
nsa=answer(10);
nsb=answer(11);

nbah=alphah*(1-iah*Xa(omegah));
nbal=alphal*(1-ial*Xa(omegal));
nba=nbah+nbal;
nbbh=alphah*(1-ibh*Xb(omegah));
nbbl=alphal*(1-ibl*Xb(omegal));
nbb=nbbh+nbbl;

sigmaah=nbah/nba;
sigmaal=nbal/nba;
sigmabh=nbbh/nbb;
sigmabl=nbbl/nbb;

mba=phi*(nsa/nba)^nu;
msa=phi*(nba/nsa)^(1-nu);
mbb=phi*(nsb/nbb)^nu;
msb=phi*(nbb/nsb)^(1-nu);

auxh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
auxl=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
gh=1+auxh*gamma_a*(iah*Xa(omegah)+(1-iah*Xa(omegah))*mba)+auxh*gamma_b*(ibh*Xb(omegah)+(1-ibh*Xb(omegah))*mbb);
gl=1+auxl*gamma_a*(ial*Xa(omegal)+(1-ial*Xa(omegal))*mba)+auxl*gamma_b*(ibl*Xb(omegal)+(1-ibl*Xb(omegal))*mbb);

g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
rt=1/rm-1+delta;
w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));

oa=(msa*theta*(sigmaah*v1h+sigmaal*v1l)*gamma_a)/(1-(1-msa*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));
ob=(msb*theta*(sigmabh*v1h+sigmabl*v1l)*gamma_b)/(1-(1-msb*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));

zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
zhbar=a*zlbar;    

rd_exp_ah=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(iah^(1+rho_a))/(1+rho_a);
rd_exp_al=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(ial^(1+rho_a))/(1+rho_a);

rd_exp_bh=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibh^(1+rho_b))/(1+rho_b);
rd_exp_bl=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibl^(1+rho_b))/(1+rho_b);

rd_exp_h=rd_exp_ah+rd_exp_bh;
rd_exp_l=rd_exp_al+rd_exp_bl;

basic_h_share=rd_exp_bh/rd_exp_h;
basic_l_share=rd_exp_bl/rd_exp_l;

basic_share=(alphah*rd_exp_bh+alphal*rd_exp_bl)/(alphah*rd_exp_h+alphal*rd_exp_l);

zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));

Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));

rd_Y_h=rd_exp_h/Yh;
rd_Y_l=rd_exp_l/Yl;

patent_count=alphah*(iah+ibh)+alphal*(ial+ibl);
trans_share=(alphah*iah*(1-Xa(omegah))+alphal*ial*(1-Xa(omegal)))*(msa*(1-((1-msa)*delta_p)^10))/((1-(1-msa)*delta_p)*patent_count)+(alphah*ibh*(1-Xb(omegah))+alphal*ibl*(1-Xb(omegal)))*(msb*(1-((1-msb)*delta_p)^10))/((1-(1-msb)*delta_p)*patent_count);


if v1h*gamma_a >= oa && v1l*gamma_a>=oa
    disp('Both types keep applied patents');
elseif v1h*gamma_a > oa && v1l*gamma_a<oa
    disp('High type keeps applied patents; low type sells applied patents');
else
    disp('Both types sell applied patents');
end

if v1h*gamma_b >= ob && v1l*gamma_b>=ob
    disp('Both types keep basic patents');
elseif v1h*gamma_b > ob && v1l*gamma_b<ob
    disp('High type keeps basic patents; low type sells basic patents');
else
    disp('Both types sell basic patents');
end

%--------------------------------------------------------------------------
%                            Generate Tables                              %
%--------------------------------------------------------------------------

phi_initial=phi; theta_initial=theta; mu_initial=mu; iota_initial=iota; gamma_a_initial=gamma_a; gamma_b_initial=gamma_b; rho_a_initial=rho_a; rho_b_initial=rho_b; chi_b_initial=chi_b; ss_initial=ss;
omegah_initial=omegah; omegal_initial=omegal; omega_avg_initial=alphah*omegah+alphal*omegal; rd_Y_h_initial=rd_Y_h; rd_Y_l_initial=rd_Y_l; trans_share_initial=trans_share; g_avg_initial=g_avg; basic_share_initial=basic_share;

% Create Appendix Table 18: Model Fit for Key Moments in the Initial Balanced Growth Path

Targets = {
    'Economic Growth Rate (1981–1985)';
    'R&D Cost/Sales of H Firms (1981–1985)';
    'R&D Cost/Sales of L Firms (1981–1985)';
    'Basic R Share of H Firms (1981–1985)';
    'Basic R Share of L Firms (1981–1985)';
    'Avg. Number of Industries of H Firms (1981–1985)';
    'Avg. Number of Industries of L Firms (1981–1985)';
    'The Share of Patents Transacted (1983)'
};

Data = {g_target^(zeta/(zeta+lambda))-1; rd_sales_h; rd_sales_l; basic_h; basic_l;  range_h; range_l; trans_target};

Model = {g_avg^(zeta/(zeta+lambda))-1; rd_Y_h; rd_Y_l; basic_h_share; basic_l_share; omegah; omegal; trans_share};

Table18 = table(Targets, Data, Model, ...
    'VariableNames', {'Targets', 'Data', 'Model'});
disp(Table18);
writetable(Table18, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table18.xlsx");

%--------------------------------------------------------------------------
%                     Calibrating to the Ending BGP                      %
%--------------------------------------------------------------------------

ss=0.24; % R&D Tax Credit in 1996-2000

% Calibrate key parameters

par_vec=@(vec)par_solver_ext2(vec,g_target2,rd_sales_h2,rd_sales_l2,range_h2,range_l2,trans_target2,basic_h2,basic_l2);
vec_guess=[2 ,   0.00015  , 1  , 5 ,   0.6,  0.2 ,  0.06  ,  0.2 ];
lb=[0,0,0,0,0,0];
ub=[1E6,1E6,1E6,1E6,1E6,1];
[par_opt,V_opt,exitflag]=fmincon(par_vec,vec_guess,[],[],[],[],[],[],[],options2);

if exitflag <=0
    disp('Parameters Not Solved')
end

iota=par_opt(1);
mu=par_opt(2);
gamma_a=par_opt(3);
chi_b=par_opt(4);
rho_a=par_opt(5);
rho_b=par_opt(6);
phi=par_opt(7);
theta=par_opt(8);

gamma_b=1.6*gamma_a;

% Compute Moments in the ending BGP

guess=[0.6,0.1,0.02,0.001,11,2,1,6,1,2,0.01];
geq_vec=@(input)geq_solver_ext(input,iota,mu,gamma_a,gamma_b,chi_b,rho_a,rho_b,phi,theta);
[answer,value,exitflag] = fsolve(geq_vec, guess, options);
if exitflag ~= 1
    disp('First-Order Optimality is not Achieved')
end


iah=answer(1);
ial=answer(2);
ibh=answer(3);
ibl=answer(4);
omegah=answer(5);
omegal=answer(6);
a=answer(7);
v1h=answer(8);
v1l=answer(9);
nsa=answer(10);
nsb=answer(11);

nbah=alphah*(1-iah*Xa(omegah));
nbal=alphal*(1-ial*Xa(omegal));
nba=nbah+nbal;
nbbh=alphah*(1-ibh*Xb(omegah));
nbbl=alphal*(1-ibl*Xb(omegal));
nbb=nbbh+nbbl;

sigmaah=nbah/nba;
sigmaal=nbal/nba;
sigmabh=nbbh/nbb;
sigmabl=nbbl/nbb;

mba=phi*(nsa/nba)^nu;
msa=phi*(nba/nsa)^(1-nu);
mbb=phi*(nsb/nbb)^nu;
msb=phi*(nbb/nsb)^(1-nu);

auxh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
auxl=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
gh=1+auxh*gamma_a*(iah*Xa(omegah)+(1-iah*Xa(omegah))*mba)+auxh*gamma_b*(ibh*Xb(omegah)+(1-ibh*Xb(omegah))*mbb);
gl=1+auxl*gamma_a*(ial*Xa(omegal)+(1-ial*Xa(omegal))*mba)+auxl*gamma_b*(ibl*Xb(omegal)+(1-ibl*Xb(omegal))*mbb);

g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
rt=1/rm-1+delta;
w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));

oa=(msa*theta*(sigmaah*v1h+sigmaal*v1l)*gamma_a)/(1-(1-msa*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));
ob=(msb*theta*(sigmabh*v1h+sigmabl*v1l)*gamma_b)/(1-(1-msb*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));

zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
zhbar=a*zlbar;    

rd_exp_ah=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(iah^(1+rho_a))/(1+rho_a);
rd_exp_al=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(ial^(1+rho_a))/(1+rho_a);

rd_exp_bh=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibh^(1+rho_b))/(1+rho_b);
rd_exp_bl=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibl^(1+rho_b))/(1+rho_b);

rd_exp_h=rd_exp_ah+rd_exp_bh;
rd_exp_l=rd_exp_al+rd_exp_bl;

basic_h_share=rd_exp_bh/rd_exp_h;
basic_l_share=rd_exp_bl/rd_exp_l;

basic_share=(alphah*rd_exp_bh+alphal*rd_exp_bl)/(alphah*rd_exp_h+alphal*rd_exp_l);

zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));

Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));

rd_Y_h=rd_exp_h/Yh;
rd_Y_l=rd_exp_l/Yl;

patent_count=alphah*(iah+ibh)+alphal*(ial+ibl);
trans_share=(alphah*iah*(1-Xa(omegah))+alphal*ial*(1-Xa(omegal)))*(msa*(1-((1-msa)*delta_p)^10))/((1-(1-msa)*delta_p)*patent_count)+(alphah*ibh*(1-Xb(omegah))+alphal*ibl*(1-Xb(omegal)))*(msb*(1-((1-msb)*delta_p)^10))/((1-(1-msb)*delta_p)*patent_count);

if v1h*gamma_a >= oa && v1l*gamma_a>=oa
    disp('Both types keep applied patents');
elseif v1h*gamma_a > oa && v1l*gamma_a<oa
    disp('High type keeps applied patents; low type sells applied patents');
else
    disp('Both types sell applied patents');
end

if v1h*gamma_b >= ob && v1l*gamma_b>=ob
    disp('Both types keep basic patents');
elseif v1h*gamma_b > ob && v1l*gamma_b<ob
    disp('High type keeps basic patents; low type sells basic patents');
else
    disp('Both types sell basic patents');
end

%--------------------------------------------------------------------------
%                            Generate Tables                              %
%--------------------------------------------------------------------------

phi_ending=phi; theta_ending=theta; mu_ending=mu; iota_ending=iota; gamma_a_ending=gamma_a; gamma_b_ending=gamma_b; rho_a_ending=rho_a; rho_b_ending=rho_b; chi_b_ending=chi_b; ss_ending=ss;
omegah_ending=omegah; omegal_ending=omegal; omega_avg_ending=alphah*omegah+alphal*omegal; rd_Y_h_ending=rd_Y_h; rd_Y_l_ending=rd_Y_l; trans_share_ending=trans_share; g_avg_ending=g_avg; basic_share_ending=basic_share;


% Create Appendix Table 19: Model Fit for Key Moments in the Ending Balanced Growth Path

Targets = {
    'Economic Growth Rate (1996–2000)';
    'R&D Cost/Sales of H Firms (1996–2000)';
    'R&D Cost/Sales of L Firms (1996–2000)';
    'Basic R Share of H Firms (1996–2000)';
    'Basic R Share of L Firms (1996–2000)';
    'Avg. Number of Industries of H Firms (1996–2000)';
    'Avg. Number of Industries of L Firms (1996–2000)';
    'The Share of Patents Transacted (2000)'
};


Data = {g_target2^(zeta/(zeta+lambda))-1; rd_sales_h2; rd_sales_l2; basic_h2; basic_l2;  range_h2; range_l2; trans_target2};

Model = {g_avg^(zeta/(zeta+lambda))-1; rd_Y_h; rd_Y_l; basic_h_share; basic_l_share; omegah; omegal; trans_share};

Table19 = table(Targets, Data, Model, ...
    'VariableNames', {'Targets', 'Data', 'Model'});
disp(Table19);
writetable(Table19, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table19.xlsx");

% Create Appendix Table 17: Parameter Values of the Extended Model
Parameter = {'γb/γa';'χa';'Xa(ω)';'Xb(ω)';'γa';'χb';'1+ρa';'1+ρb';'γa';'χb';'1+ρa';'1+ρb'};

Description = {
    'Step Size Ratio';
    'Applied R Cost, Scale';
    'Applied R, Within-Scope Prob.';
    'Basic R, Within-Scope Prob.';
    'Applied R Step Size (Initial BGP)';
    'Basic R Cost, Scale (Initial BGP)';
    'Applied R Cost, Elasticity (Initial BGP)';
    'Basic R Cost, Elasticity (Initial BGP)';
    'Applied R Step Size (Ending BGP)';
    'Basic R Cost, Scale (Ending BGP)';
    'Applied R Cost, Elasticity (Ending BGP)';
    'Basic R Cost, Elasticity (Ending BGP)'
};

Value = {
    gamma_b/gamma_a;                           % γb/γa
    chi_a;                                     % χa
    "exp(-3.837) * |omega|^0.602";             % Xa(ω)
    "exp(-4.944) * |omega|^0.932";             % Xb(ω)
    gamma_a_initial;                                   % γa
    chi_b_initial;                                     % χb
    1+rho_a_initial;                                   % 1+ρa
    1+rho_b_initial;                                    % 1+ρb
    gamma_a_ending;                                   % γa
    chi_b_ending;                                     % χb
    1+rho_a_ending;                                   % 1+ρa
    1+rho_b_ending                                    % 1+ρb
};

Identification = {
    'Akcigit et al. (2021)';
    'Normalization';
    'Regression';
    'Regression';
    'Growth Rate';
    'Basic Research Share';
    'R&D Cost/Sales';
    'Ratio (H and L)';
    'Growth Rate';
    'Basic Research Share';
    'R&D Cost/Sales';
    'Ratio (H and L)'
};


Table17 = table(Parameter, Description, Value, Identification, ...
    'VariableNames', {'Parameter', 'Description', 'Value', 'Identification'});

disp(Table17);
writetable(Table17, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table17.xlsx");


% Create Table 10: Effects of Key Parameters

Decomposition_extension